%--------------------------------------------------------------------------
%----- 2-D DIFFUSION-BASED RESPONSE MATRIX SOLVER
%----- TWO-DIMENSIONAL BENCHMARK PROBLEMS:
%----- 1. IAEA 
%----- 2. BIBLIS
%----- 3. KOEBERG (4 group)

name = 'two_d_benchmarks';
%    0.019487941261977
% 
% max(abs(rel))
% 
% ans =
% 
%    0.057260021353619

prob      = 0;  % 1, 2, or 3 for IAEA BIBLIS or KOEBERG
locsolv   = 0;  % 0 for KERNEL or 1 for FINITE DIFFERENCE
locfit    = 0;
reflector = 0;
order = 0; % legendre order
bc(1) = 1; % (L)
bc(2) = 1; % (R)
bc(3) = 1; % (B)
bc(4) = 1; % (T)
abdat=0;
if prob==0
    
%--- MATERIALS
numtypes    = 4;   % number of element types
numg        = 2;   % number of groups
numm        = 4;   % number of materials
dat = [  1.4493    0.0250    0.0079    1.0000    0.0000    0.0000 % fresh
         0.3807    0.1042    0.1692    0.0000    0.0151    0.0000 
         1.4479    0.0258    0.0069    1.0000    0.0000    0.0000 % once
         0.3708    0.1200    0.1745    0.0000    0.0148    0.0000
         1.4494    0.0262    0.0060    1.0000    0.0000    0.0000 % twice
         0.3676    0.1191    0.1625    0.0000    0.0147    0.0000 
         1.3200    0.0257622 0.0000    1.0000    0.0000    0.0000 % ref (biblis)
         0.2772    0.0715960 0.0000    0.0000    0.0231060 0.000];
  
geo = 1;

if geo == 1
    
    elements = [ 2 2 2 
                 2 1 2
                 2 2 2];
   
elseif geo == 0
     
elements = [    3 1 2 1 2 1 2 3 0
                1 2 1 2 1 2 1 3 0
                2 1 2 1 2 1 3 3 0
                1 2 1 2 1 2 3 0 0
                2 1 2 1 2 3 3 0 0
                1 2 1 2 3 3 0 0 0 
                2 1 3 3 3 0 0 0 0
                3 3 3 0 0 0 0 0 0 
                0 0 0 0 0 0 0 0 0               
          ];      
elements(:, :) =1;
elements(elements==0)=0;
      
elseif geo == 2      
elements = [   
                3 1 1 2 2 1 1 2 2 1 1 2 2 3 3 4 4
                1 2 2 1 1 2 2 1 1 2 2 1 1 3 3 4 4
                1 2 2 1 1 2 2 1 1 2 2 1 1 3 3 4 4
                2 1 1 2 2 1 1 2 2 1 1 3 3 3 3 4 4
                2 1 1 2 2 1 1 2 2 1 1 3 3 3 3 4 4
                1 2 2 1 1 2 2 1 1 2 2 3 3 4 4 4 4
                1 2 2 1 1 2 2 1 1 2 2 3 3 4 4 4 4
                2 1 1 2 2 1 1 2 2 3 3 3 3 4 4 4 4
                2 1 1 2 2 1 1 2 2 3 3 3 3 4 4 4 4
                1 2 2 1 1 2 2 3 3 3 3 4 4 4 4 4 4 
                1 2 2 1 1 2 2 3 3 3 3 4 4 4 4 4 4 
                2 1 1 3 3 3 3 3 3 4 4 4 4 4 4 4 4
                2 1 1 3 3 3 3 3 3 4 4 4 4 4 4 4 4
                3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
                3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
                4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4               
                4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4   ]; 
%elements(:, :) =1;
elements(elements==4)=0;

elseif geo==3           
elements =[ 3 3 3 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4
            3 3 3 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4
            3 3 3 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 3 3 3 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 3 3 3 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 3 3 3 3 3 3 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 4 4 4 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4
            2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4
            1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 
            1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 
            1 1 1 2 2 2 1 1 1 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 
            2 2 2 1 1 1 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4
            2 2 2 1 1 1 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4
            2 2 2 1 1 1 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ];
elseif geo==4      
elements = [3 3 3 3 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
            3 3 3 3 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
            3 3 3 3 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
            3 3 3 3 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 4 4 4 4            
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
            1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 
            2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            2 2 2 2 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4];
end
      %elements=elements*0.0 + 1;
xcm  = [   0     21  ]/geo;
xfm  = [     1       ]*floor(10/geo);  
ycm  = [   0     21  ]/geo;
yfm  = [     1       ]*floor(10/geo);  
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ];     
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ]; 
abdat=0; % not used  



elseif prob==1
%--------------------------------------------------------------------------
% PROBLEM 1: IAEA
%--------------------------------------------------------------------------
%--- MATERIALS
numtypes    = 4;   % number of element types
numg        = 2;   % number of groups
numm        = 4;   % number of materials
buck = 0.8e-4;
dat = [ 1.500   0.030   0.000  1.000  0.000  0.000  
        0.400   0.080   0.135  0.000  0.020  0.000
        1.500   0.030   0.000  1.000  0.000  0.000
        0.400   0.085   0.135  0.000  0.020  0.000
        1.500   0.030   0.000  1.000  0.000  0.000  
        0.400   0.130   0.135  0.000  0.020  0.000 
        2.000   0.040   0.000  0.000  0.000  0.000
        0.300   0.010   0.000  0.000  0.040  0.000 ];       
abdat=0; % not used  

%--- PLACE THE ELEMENTS            
elements = [3 2 2 2 3 2 2 1 0
            2 2 2 2 2 2 2 1 0
            2 2 2 2 2 2 1 1 0
            2 2 2 2 2 2 1 0 0
            3 2 2 2 3 1 1 0 0
            2 2 2 2 1 1 0 0 0 
            2 2 1 1 1 0 0 0 0
            1 1 1 0 0 0 0 0 0 
            0 0 0 0 0 0 0 0 0
          ];
% elements = [3 2 2 2 3 2 2 1 3
%             2 2 2 2 2 2 2 1 3
%             2 2 2 2 2 2 1 1 3
%             2 2 2 2 2 2 1 3 3
%             3 2 2 2 3 1 1 3 3
%             2 2 2 2 1 1 3 3 3 
%             2 2 1 1 1 3 3 3 3
%             1 1 1 3 3 3 3 3 3 
%             3 3 3 3 3 3 3 3 3
%           ];      
      %elements=elements*0.0 + 1;
xcm  = [   0     20  ];
xfm  = [     1       ]*10;  
ycm  = [   0     20  ];
yfm  = [     1       ]*10;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ]; 

elseif prob ==2
    
%--------------------------------------------------------------------------
% PROBLEM 2: Biblis
%--------------------------------------------------------------------------
%--- MATERIALS
numtypes    = 8;   % number of element types
numg        = 2;   % number of groups
numm        = 8;   % number of materials
dat = [ 1.4360000  0.0272582  0.0058708  1.0000000  0.0000000  0.0000000
        0.3635000  0.0750580  0.0960670  0.0000000  0.0177540  0.0000000
        1.4366000  0.0272995  0.0061908  1.0000000  0.0000000  0.0000000
        0.3636000  0.0784360  0.1035800  0.0000000  0.0176210  0.0000000
        1.3200000  0.0257622  0.0000000  0.0000000  0.0000000  0.0000000
        0.2772000  0.0715960  0.0000000  0.0000000  0.0231060  0.0000000
        1.4389000  0.0274640  0.0074527  1.0000000  0.0000000  0.0000000
        0.3638000  0.0914080  0.1323600  0.0000000  0.0171010  0.0000000
        1.4381000  0.0272930  0.0061908  1.0000000  0.0000000  0.0000000
        0.3665000  0.0848280  0.1035800  0.0000000  0.0172900  0.0000000
        1.4385000  0.0273240  0.0064285  1.0000000  0.0000000  0.0000000
        0.3665000  0.0873140  0.1091100  0.0000000  0.0171920  0.0000000
        1.4389000  0.0272900  0.0061908  1.0000000  0.0000000  0.0000000
        0.3679000  0.0880240  0.1035800  0.0000000  0.0171250  0.0000000
        1.4393000  0.0273210  0.0064285  1.0000000  0.0000000  0.0000000
        0.3680000  0.0905100  0.1091100  0.0000000  0.0170270  0.0000000];
abdat=0; % not used  
%--- PLACE THE ELEMENTS   
elements = [ 1 8 8 2 2 6 6 1 1 7 7 1 1 4 4 3 3   % 1
             8 1 1 8 8 2 2 8 8 1 1 1 1 4 4 3 3   % 2
             8 1 1 8 8 2 2 8 8 1 1 1 1 4 4 3 3   % 3
             2 8 8 1 1 8 8 2 2 7 7 1 1 4 4 3 3   % 4
             2 8 8 1 1 8 8 2 2 7 7 1 1 4 4 3 3   % 5
             6 2 2 8 8 2 2 8 8 1 1 8 8 4 4 3 3   % 6
             6 2 2 8 8 2 2 8 8 1 1 8 8 4 4 3 3   % 7           
             1 8 8 2 2 8 8 2 2 5 5 4 4 3 3 3 3   % 8
             1 8 8 2 2 8 8 2 2 5 5 4 4 3 3 3 3   % 9
             7 1 1 7 7 1 1 5 5 4 4 4 4 3 3 3 3   % 10
             7 1 1 7 7 1 1 5 5 4 4 4 4 3 3 3 3   % 11
             1 1 1 1 1 8 8 4 4 4 4 3 3 3 3 3 3   % 12
             1 1 1 1 1 8 8 4 4 4 4 3 3 3 3 3 3   % 13       
             4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3   % 14
             4 4 4 4 4 4 4 3 3 3 3 3 3 3 3 3 3   % 15
             3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   % 16
             3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3   % 17
             ];    
xcm  = [   0     20.1226/2  ]*2;
xfm  = [     1       ]*2;
ycm  = [   0     20.1226/2  ]*2;
yfm  = [     1       ]*100;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);
% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ];    
% ELEMENT 5 - material specification
mt(:,:,5)   = [  5 ]; 
% ELEMENT 6 - material specification
mt(:,:,6)   = [  6  ]; 
% ELEMENT 7 - material specification
mt(:,:,7)   = [  7 ]; 
% ELEMENT 8 - material specification
mt(:,:,8)   = [  8  ];    

end



tol     = [1.e-9 1.e-9]; % tolerance for norms             
maxit   = 20;

swpi = 1; % currently not used
swnt = 1; % currently not used
swgm = 0; % ''

%-------------------------------------------------------------------------
% RUN SOLVER and OUTPUT EITHER TO FILE OR SCREEN      

ord  = [2 1 0];
% mesh = [150 125 100 75 50 25 20 15 10 5 4 3 2];

%for i = 1:length(ord)
% 
%order = ord(i);
input = makeinput( name, numtypes, numg, numm, order, dat, abdat, ...
                   elements, bc, xcm, xfm, ycm, yfm, mt, tol, maxit, ... 
                   swpi, swnt, swgm, 1 ) ; 
%input.x0 = x1;
input.locsolv   = locsolv;
input.locfit    = locfit;
input.reflector = reflector;
%input.x0(1:end-2) = 1.1769413;
tic
[x1,normk1,it,nrm1,fr] = solverPI_mc(input);

xx = x1(1:end-2);
xx = reshape(xx, 8, length(xx)/8)

disp(['  iteration number: ', num2str(it), ' norm(x) = ', num2str(nrm1(it))])
fprintf(1,'  final keff = %12.8f\n',x1(end-1));
toc

%  FR4(:,i) = fr;
%  k4(i)    = x1(end-1);
 
% end


if locsolv==0 && locfit==1
    frf = fr;
    rel2=(frf-frr)./frr; rel2(isnan(rel2))=0;
    kf  = x1(end-1);
elseif locsolv==0 && locfit==0
    frm = fr;
    rel1=(frm-frr)./frr; rel1(isnan(rel1))=0;
    km  = x1(end-1);
else
    frr = fr;
    kr  = x1(end-1);
end


